package scu.maqiang.piezoelectricity;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T3PE;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.T3Type;
import scu.maqiang.numeric.DirectSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

import java.util.Arrays;

public class AxisPiezoelectricity2DT3 {

	public static void main(String[] args) {
		long begin = System.nanoTime();
		Mesh2T3 mesh = new Mesh2T3();
		mesh.square2D(100, 100, T3Type.Left);
		mesh.translate(1.0, 0.0);
		FES2T3PE fs = new FES2T3PE(mesh);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		double[] RHS = new double[fs.GetNdof()];
		double[] x = RHS.clone();
		double[][] AmatD = new double[][] {{14.9e10, 9.8e10,  9.8e10,    0.0,     0.0,     -5.4}, 
			                              {9.8e10,  14.3e10,  9.8e10,    0.0,     0.0,     -5.4},
			                              {9.8e10,   9.8e10, 14.3e10,    0.0,     0.0,     -5.4},
			                              {0.0,         0.0,     0.0, 2.2e10,    11.7,      0.0}, 
			                              {0.0,         0.0,     0.0,  -11.7,  8.0e-9,      0.0},
			                              {5.4,         5.4,     5.4,    0.0,     0.0,   7.2e-9}};
	    fs.assembleStiff(new double[][][] {AmatD}, BVPType.CONSTITUTE_MATRIX_AXIS, A);
		//System.out.println(A);
	    fs.assemblePESource(new double[] {100.0, 0.0, 200}, BVPType.AXISYMMETRIC, RHS);
	    fs.applyBC_MBN(A, RHS, new int[] {0, 1}, new double[] {0.0, 0.0}, 1, 2, 3, 4);
	    fs.applyBC_MBN(A, RHS, new int[] {2}, new double[] {0.0}, 1, 2, 3);
	    
	    DirectSSolver dsolver = new DirectSSolver(A);
		double[][] uvphi = new double[3][mesh.getNv()];
	    x = dsolver.LUSolve(RHS, 0, 1e-16);
		fs.extractComponent(x, uvphi);
		System.out.println("SOlution: ");
		System.out.println("u min:" + MVO.min(uvphi[0]) + "\tmax: " + MVO.max(uvphi[0]));
		System.out.println("v min:" + MVO.min(uvphi[1]) + "\tmax: " + MVO.max(uvphi[1]));
		System.out.println("phi min:" + MVO.min(uvphi[2]) + "\tmax: " + MVO.max(uvphi[2]));
	    //A.save("SRMatrixFile.txt");
	    //MVO.saveArray("RHSFile.txt", RHS);
		Arrays.fill(x, 0.0);
	    NewIterSSolver solver = new NewIterSSolver(A);
		solver.setEps(1.0e-16);
		solver.PCGSSSOR(RHS, x, 1.5, 1);
	    //solver.setEps(1.0e-10);
	    //solver.GMRES_SSOR(RHS, x, 1.5, 200, 1);
	    fs.extractComponent(x, uvphi);
	    mesh.toTecplot("PESolution2D.dat", uvphi);
	    System.out.println("SOlution: ");
	    System.out.println("u min:" + MVO.min(uvphi[0]) + "\tmax: " + MVO.max(uvphi[0]));
	    System.out.println("v min:" + MVO.min(uvphi[1]) + "\tmax: " + MVO.max(uvphi[1]));
	    System.out.println("phi min:" + MVO.min(uvphi[2]) + "\tmax: " + MVO.max(uvphi[2]));
		long end = System.nanoTime();
		System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
	}

}
